On the reliability of recent Monte Carlo studies of dilute systems of localized spins 

interacting with itinerant carriers. 
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In this paper, we discuss magnetic properties of dilute systems of localized spins interacting with 
itinerant carriers. More precisely, we compare recent available Monte Carlo results with our two 
step approach (TSA) calculations. The TSA consists first on the determination of the magnetic 
couplings and then on a proper treatment of the resulting effective dilute Heisenberg Hamiltonian. 
We show important disagreement between the Monte Carlo results (in principle exact) and our TSA 
calculations. We analyze the origins and shed light on the reasons of those dissensions. In contrast to 
one could expect, we demonstrate that the available MC calculations suffer from severe numerical 
shortcomings. More precisely, (i) the finite size effects appear to be huge in dilute systems, (ii) 
the statistical sampling (disorder configurations) was far too small, and (iii) the determination of 
the Curie temperature was too rough. In addition, we provide new arguments to explain a recent 
disagreement between the Monte Carlo simulations and TSA for the model study of the well known 
III-V diluted magnetic semiconductor Gai^^Mn^As. We hope that this work will motivate new 
systematic large scale Monte Carlo calculations. 
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I. INTRODUCTION 

During the last decade, the study of the effects of 
disorder on magnetic properties (Curie/Neel tempera- 
tures, magnetic excitations spectrum, spin susceptibil- 
ity,...) of disordered/dilute magnetic materials became a 
central field of research for both theoreticians and ex- 
perimentalists. It is now clear that even for a qual- 
itative understanding of the underlying physics, disor- 
der, dilution and thermal fluctuations should be handled 
in a proper way. Among widely studied materials one 
finds oxides manganite as for example Lai-^Sr^MnOa 
or Lai-^Ba^MnOs, which exhibit the fascinating giant 
magneto-resistance phenomenon^! 3 -! 4 -. The diluted mag- 
netic semiconductors (DMS) as the well known III-V 
compound Ga, 1 - x Mn x As£&Z& or the II- VI Zni^Cr^Te^ 
have attracted a considerable attention. The inter- 
est for DMS is partly motivated by their technolog- 
ical potential for spintronics. One can also mention 
the so-called d° materials as the wide band gap ox- 
ides HfC<2, ZrC>2, CaO, ZnO, or the irradiated graphite 
layers,.. i 10 ' n i 12 ' 13 i 14 . These new materials may become 
an alternative for spintronic applications. The Hcusler 
alloys as Ni2+ a Mni_ a; S n 15 i :L6 or the double perovskites 
as S^FeMoO o 17 ' 18 are also materials in which disorder 
clearly plays a crucial role. 

Ab-initio based studies have provided the most effi- 
cient and reliable tool, to allow theoretical and quan- 
titative study of the magnetic properties without ad- 
justing parameters. As an example, it was possible to 
study in great details and quantitatively the magnetic 



properties of the III-V diluted magnetic semiconductor 
Gai-^Mn^As, both in the presence or absence of addi- 
tional compensating defect a 19 ' 20 i 21 . In order to test the 
interest of a particular material for technological appli- 
cations, the ab initio approach provides the most appro- 
priate tool. However, because it is material specific and 
very complex, it does not allow to draw general conclu- 
sions which could be valid for a whole family of mate- 
rials. The study of relevant minimal model is crucial 
and allows to fill this gap. The model approach is essen- 
tial to understand the influence of a particular physical 
free parameter. It can also help for the determination of 
potential spintronic candidates. The minimal one band 
model that describes itinerant carriers (holes/electrons) 
interacting with localized moments reads, 



^piJSi ■ Si + 2J1 



(1) 



In the first term, ty=t if i and j are nearest neighbor, 
otherwise it is 0. The random variable pi is one if the 
site is occupied by a magnetic impurity, otherwise it is 0, 
we denote x the concentration of impurities. Sj denotes 
the localized impurity spin at site i. For most of the 
materials S is large enough, we consider here the case of 
classical spins only. Moreover in ferromagnets the quan- 
tum fluctuations are not relevant. In the case of man- 
ganites pi = 1 (one Mn per site) Ji is the Hund coupling 
(Jh ) between the itinerant e g carriers and the localized 
tig magnetic moments (S = 3/2). In diluted magnetic 
semiconductors, Ji is the local p-d coupling (J p( j). The 
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additional last term describes the effects of the disorder 
induced by the substitution of a cation by another one. 
Note that this term includes in particular the electro- 
static contribution which originates from the difference 
of charge between the substituted cation and the original 
one. 

In this manuscript, all the calculations are performed 
on a simple cubic lattice. We also neglect in eq.JT]) the 
last term. The case of finite on-site potential will be dis- 
cussed elsewhere^. We expect the more realistic case of 
a multi-band Hamiltonian to lead to similar conclusions 
to those that will be presented in this manuscript. 



II. THE TWO STEPS APPROACH 
PREDICTIONS 

Reliable calculations requires several crucial conditions 
to be fulfilled. In other words, one should treat (i) the 
local coupling in a non perturbative way, (ii) the ther- 
mal and transverse spin fluctuations properly, and (iii) 
the disorder and dilution effects should be treated ex- 
actly (beyond an effective medium approach). The last 
condition allows for the localization of the itinerant car- 
riers. This will strongly affect both magnetic and trans- 
port properties. The full Monte Carlo method is a pri- 
ori the best procedure, since localized spins degrees of 
freedom and itinerant carriers are treated exactly, simul- 
taneously and on equal footings. Note that the word 
full is added to avoid confusion with the Monte Carlo 
treatment of the effective Heisenberg Hamiltonian in the 
two step approach procedure presented below. An alter- 
native approach which fulfill these requirements is the 
Two Step Approach (TSA). TSA is in spirit similar to 
what is widely used in the ab-initio based studies. It 
consists in two different stages. First, for a given con- 
figuration of disorder, at T — K, one diagonalizes the 
Hamiltonian given by Eq. (TJJ.This provides the spin re- 
solved one particle Green's Function of the itinerant car- 
rier (GF) Gfj(u>) (where a =f or J.). In this first step the 
appropriate magnetic texture for the localized spins is 
used. For the fcrromagnctism. all localized spins arc par- 
allel. Note that, the calculations are done without any 
approximations at T = K , thus vertex corrections and 
the multiple scattering of the carriers on the magnetic 
impurities are included exactly. From Gfj(u>) one then 
calculate the magnetic couplings between the localized 
moments (integrate out the carriers degrees of freedom) . 
We end up with the disordered/dilute Heisenberg Hamil- 
tonian which reads, 



(2) 



where the magnetic couplings are given by the well known 

2S 24 

expressio n i 

1 f +0 ° 

J ij = -^ m (_ /(w)A i Gj J . (w) Aj-Gj^wjdw) (3) 



In the present case Ai=Aj=J p d (exchange splitting). 
The second step of the TSA consists in diagonalizing 
reliably the effective disordered Heisenberg Hamiltonian 
(Eq. ©). Note that the diagonalization can be per- 
formed either by a Monte Carlo treatment (as mentioned 
above) or within the semi-analytical method called Self- 
Consistent Local RPA (SC-LRPA)^. Once again, to 
avoid any confusion, we underline and distinguish be- 
tween the full MC study (no mapping to Heisenberg 
Hamiltonian) and the Monte Carlo that can be used to 
diagonalize Hamiltonian (JSJ in the TSA. It is important 
to recall that TSA have already been successfully used for 
the study of various realistic materials. In these studies, 
the first step of the TSA was provided by ab initio cal- 
culations (realistic exchange integrals). The agreement 
between calculated Curie temperatures and the experi- 
mental data was excellen t 25 ! 26 ! 27 ' 28 . Concerning the sec- 
ond step, it has also been shown that for a given set 
of exchange integrals, the SC-LRPA is as accurate as 
Monte Carlo calculation s] 29 ' 30 ' 31 . As it will be seen in the 
following, one great advantage of the TSA is the very 
large system sizes that can be handled, but are still in- 
accessible within the full Monte Carlo method. Thus, in 
this manuscript the second step of the TSA is performed 
within SC-LRPA. 
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FIG. 1: (Color online.) Nearest neighbor coupling as a func- 
tion of 1/N (where N is the number of sites) for different 7 
(the hole density is — jx) in the case of a non dilute (e.g. 
x = 1) and non disordered system (V = 0), in the double 
exchange limit (JS=oo). 



A. Results for non dilute systems 

Before discussing the case of dilute magnetic systems, 
it is important to discuss the non dilute case, e.g x = 1. A 
comparison between the full Monte Carlo treatment and 
the two steps approach has been done in the limit JS= 00 
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(double exchange (DE) limit )2S. In this case, the near- 
est neighbor coupling between localized moments reads 
JijS 2 = Uj(c\cj), the spin index a is irrelevant in this 
case, the fermion operators correspond to the majority 
spin. The comparison was performed in the presence of 
Anderson disorder, e.g the variables ej in Eq.([T]) were cho- 
sen randomly in the interval [— V/2,V/2], V denotes the 
strength of the disorder. Even in the presence of disor- 
der an excellent agreement between the full Monte Carlo 
calculations and the two step approach was achieved, the 
difference in the Curie temperature was within 10 %. We 



ligiblc. Similar conclusions have been reached in the case 
where the on-site disorder (V ^ in the previous Hamil- 
tonian)is also included^. We have even obtained that the 
average over only few configurations of disorder is enough 
to provide an accurate value. Thus, one can understand 
the success of the full Monte Carlo simulations to provide 
accurate Curie temperature for the study of the non di- 
lute DE model even though in these studies the systems 
sizes are typically of the order of N w 5 3 sites. The re- 
striction to such small clusters is in fact due to the very 
large CPU time and memory costs of the standard full 
Monte Carlo simulations. 
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FIG. 2: (Color online.) Curie temperature (SC-LRPA) ob- 
tained within the two step approach as a function of 7 for dif- 
ferent system size, in the case of a non dilute system (x = 1) 
and in the double exchange regime (JS=oo) and V = 0. 

now show that finite size effects are very small for the 
non dilute case. This will explain the agreement found 
between TSA calculations performed on large systems 
and full Monte Carlo simulations done with much smaller 
clusters. In Fig. [TJ the nearest neighbor magnetic cou- 
pling (Ji) between impurities as a function of 1/N (N is 
the number of sites) is plotted. It is clearly observed that 
Ji depends very weakly on the system size. The finite size 
effects are negligible at any carrier density (p^ = jx). 
In Fig. [2J we have plotted the Curie temperature as a 
function of the carrier density per magnetic moment (7). 
We recall that in the clean limit the local SC-RPA re- 
duces to the standard RPA calculation. In the case of 
classical spins the Curie temperature is given by Tc = 
f (~k Eq^o the magnon excitation dispersion 

reads E(q) = JiS 2 (l — ^(cos(q x a)+cos(q y a)+cos(q z a))). 
We clearly observe that the finite size effects are com- 
pletely negligible even for very small systems. For the 
smallest system (N = 4 3 ), the Curie temperature is al- 
ready accurate within less than 2%. Thus, in the case 
of clean systems and in the double exchange regime 
(JS=oo), the finite size effects arc indeed completely neg- 



FIG. 3: (Color online.) Spin resolved Density of states for 
a fixed concentration of magnetic impurities x = 0.065 and 
different values of the local coupling JS. The spin down (resp. 
up) band is represented by the continuous (resp. dotted) line. 
Energies are in unit of the bandwidth W = 12t. 

However, as it will be seen in the following subsection, 
for dilute systems, the use of too small system sizes as 
well as too weak statistical sampling affect the results in 
a dramatic way. 

B. Results for dilute systems 

For simplicity and to allow a direct comparison with 
the available full Monte Carlo data, we set the impurity 
concentration to x — 0.065. We will vary both, the am- 
plitude of J (its sign is irrelevant for classical spins) and 
the carrier density. Fig.[3]shows the DOS for various cou- 
pling JS. We observe that the impurity band (IB) splits 
from the valence band for values of JS > W. In the 
strong coupling regime, JS = 2W, the IB is completely 
separated from the valence band, the magnetic couplings 
are expected to be dominated by the double exchange 
mechanism. Fig. [¥] shows the couplings as a function of 
the impurity distance, for two different carrier concen- 
trations and different values of JS. Note that, the cal- 
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FIG. 4: (Color online.) Magnetic couplings (in unit of t) as 
a function of the distance between magnetic impurities. The 
density of impurity is x — 0.065 and the hole concentration 
is respectively (a) pn = 0.125a; (or 7 = 0.125) and (b) pn = 
0.25a: (or 7 = 0.25). The magnetic exchanges are calculated 
for different values of JS. In both figures, the inset focuses on 
the most "relevant" couplings. 



culations were performed on a sufficiently large cluster 
(N = 20 3 ) to neglect the finite size effects. The discus- 
sion of the finite size effects will be done in a forthcoming 
paragraph. First, we observe in Fig. 0] that the couplings 
are all ferromagnetic and of relatively short range. For 
x = 0.065, the typical average distance between magnetic 
impurities is d ps 2.4a. It corresponds to the distance be- 
tween impurities on the ordered underlying super-lattice. 
For the smallest values of JS, the couplings are small 
at any distances, this should lead to very small Curie 
temperatures. As JS is increased, the nearest neighbor 
coupling increases strongly until it saturates whilst the 
range of the couplings is reducing (clearly seen in the 
insets). For very large JS (JS » W), the couplings are 
essentially reduced to the nearest neighbor coupling only, 
as expected in the double exchange limit. The limiting 
value is expected to be independent of the carrier density 
and corresponds roughly to the kinetic energy of a sin- 



gle hole hopping between two sites. Note that there are 
also some additional contributions coming from one hole 
localized on three sites, 2 holes on three sites and so on. 
One can already say that, in the large JS limit, the Curie 
temperature is expected to vanish since the average dis- 
tance between impurities is larger than the range of the 
couplings. In other words, the system is below the per- 
colation threshold. For JS « 2W (see DOS in fig.l), the 
couplings at distances near d are tiny, thus one expects 
very small Curie temperatures. 

In Fig. [5J we have plotted the Curie temperature as a 
function of JS. The carrier density is the same as used in 
Fig. [4] We again recall that the Monte Carlo treatment of 
this Hamiltonian would lead to similar results as it was al- 
ready shown in several publication o 19 ' 29 . For both values 
of 7, we observe that Tc increases until it reaches a max- 
imum, then it decreases until it vanishes beyond a critical 
value which depends on the magnetic impurity concen- 
tration. Note that the maximum is located at JS ~ 0.5 
for 7 = 0.125 and JS w 0.6 for 7 = 0.25. The high- 
est T c are respectively T™ ax (j = 0.125) w 3.10~ 4 W 
and Tg lax (j = 0.25) « 5.10~ 4 W. For small values of 
J (perturbative limit) one notices that within SC-LRPA 
Tc oc J 2 , as expected from Eq. ([3]) . In Fig. [6l we evaluate 
the role of the disorder and the importance of the thermal 
fluctuations. For that purpose, we compare Tc calcu- 
lated within SC-LRPA to the Mean Field Virtual Crys- 
tal Approximation (MF-VCA) value. We remind that 
the MF-VCA value is T c /F ~ ycvl = fx £ l#0 J 0l - Within 
this approximation both thermal and transverse fluctua- 
tions and disorder are treated at the lowest order. In fact 
the disorder (dilution)is neglected and just appears as a 
trivial prefactor (x) in the expression of Tc- In the limit 
of the very small J (inset of Fig. [6]), we observe a good 
agreement between SC-LRPA and MF-VCA. Within the 
perturbative limit, and for relatively weak densities of 
carriers only, one could indeed expect such a behavior 
(similar results are observed for 7 = 0.25 of Fig.[5ji). As 
we increase JS, we observe important quantitative and 



qualitative differences, as T 
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Let us now show that Tc is in fact mainly controlled 
by the typical coupling namely J(d), where d was defined 
previously. For x = 0.065 we remind that the average dis- 
tance between impurities is d = d ps 2.4a. In Fig. [7] we 
have plotted the variation of this coupling as a function 
of JS/W. We clearly see that the shape of the calculated 
curve is very similar to that obtained in Fig. [3^. We find 
that the ratio R = j^ s2 * s almost independent of JS 
and R ps 4. On the other hand, the non monotonic be- 
haviour of J(d) is qualitatively different from the mono- 
tonic increase of J\ as seen in the inset. The limiting 
value of J\ (for JS — » 00) is expected to be close to that 
of the clean system as plotted in Fig.l. For large JS, J\ 
dominates and leads to the unrealistic MF-VCA critical 
temperatures (see Fig- EJ) - 
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FIG. 5: (Color online.) Curie temperature within SC-LRPA 
(in units of the bandwidth W=12t) as a function of JS/W. 
The magnetic impurity concentration is fixed x — 0.065. The 
hole concentration per impurity is (a) 7 = 0.125 and (b) 7 = 
0.25. The insets represent the Monte Carlo (MC) calculations 
from ref.— done for a two band model (Nb = 2) in which (a) 
7 = 0.5 and (b) 7 = 0.25 



III. COMPARISON WITH AVAILABLE FULL 
MONTE CARLO CALCULATIONS 

In this section, we propose to compare our results with 
recent full Monte Carlo simulations results available in 
the literature. We will see that the quoted full Monte 
Carlo results arc questionable since suffering from se- 
vere numerical shortcomings. To our knowledge, the only 
available full Monte Carlo study of the diluted Hamilto- 
nian given by eq. |T]) was performed by Popcscu et al.— , 
no other group has performed similar studies in the dilute 
regime. In their MC simulations, Popescu et al. consider 
the case of two identical independent bands. Thus, both 
Tc and the carrier densities used in their simulations are 
simply divided by two in order to allow a direct compar- 



FIG. 6: (Color online.) MF-VCA and full SC-LRPA calcu- 
lations of the Curie temperature as a function of JS. The 
magnetic impurity concentration is x — 0.065 and density of 
hole per carrier is 7 = ph/x = 0.125. 



ison with our calculations (performed on a single band 
model). The MC data of ref.— are shown both in the 
insets of Fig. [5K and b. The MC critical temperatures 
were obtained for systems with N = 5 3 to 6 3 sites and 
the average over disorder was performed over 7 config- 
urations of disorder only. Let us notice that a system 
of size N = 5 3 (resp. 6 3 ) contains only 8 (resp. 14) 
magnetic impurities; thus for a density of carriers per 
impurity 7 = 0.25, such small systems respectively con- 
tain only 2 and 3 holes in the whole cluster. As seen 
in Fig. El several crucial differences appear immediately 
with respect to our results. First, we notice that for both 
carrier densities (7 = 0.125 and 0.25) the Monte Carlo 
simulations predict critical Tc much higher than those 
obtained within the two steps approach. Indeed, if one 
refers for example to Fig. [5^, one observes that the MC 
values are 15 to 60 times larger! In particular, the highest 
Tc arc typically at least 20 times larger than our results. 
For small values of JS, one notices that SC-LRPA and 
the Monte Carlo simulations do not predict the same be- 
havior (quadratic in J in our case). 

For large JS, whilst our calculated Tc? is strongly sup- 
pressed and vanishes for values of JS > J C S where 
J C S « 2 W, the MC value decreases much slowly and 
seems to saturate at large JS (see Fig. [Sk). For example, 
for 7=0.125 and JS=1.2 W (or p h =0.25 for the 2 bands 
!f c /N b = 3 ■ 10~ 3 W whilst we find 
7.510" 5 W, thus the ratio Tg c ' /T^ RPA = 40, 
which is huge! Additionally, there is an other unexpected 
and surprising result. The MF-VCA Curie temperature 
which should in principle be an upper bound is in fact 
smaller than the Monte Carlo value. This is in con- 
flict with the opposite and natural expectation. We re- 
call that the Monte Carlo method treats both disorder 
effects and thermal fluctuations exactly. For example, 
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FIG. 7: (Color online.) J(d)S 2 in unit of t as a function of 
the local coupling JS/W. The density of holes is fixed and 
set to ph = 0.125x. In the inset we have plotted the nearest 
neighbor coupling J\S 2 as a function of JS/W. 



from Fig. [6] and for JS=0.4 W, we find t^ f ~ vca ^ 4 • 
1CT 4 W and T£ /c = 50 10" 4 W. Additionally, it was al- 
ready shown several times in the literature, especially for 
the study of the diluted magnetic semiconductors, that 
r^MF-vGA i ar g e iy overestimates the real Curie temper- 
ature, this clearly put a question mark on the validity of 
the full MC results. 

On the basis of what was previously discussed, we now 
address the reasons which may explain the huge differ- 
ences with TSA and why the Mean Field VCA becomes a 
lower bound with respect to the full Monte Carlo results. 
To be more specific, it will be shown that the available 
Monte Carlo results suffer from several combined effects: 
(i) strong finite size effects, (ii) insufficient sampling (too 
few disorder configurations) and (iii) inaccurate and ap- 
proximative procedure for the determination of the Curie 
temperature. 



IV. ORIGINS OF THE DISSENSIONS 
BETWEEN TSA AND MONTE CARLO 
SIMULATIONS 

In this section we will discuss both the importance of 
the statistical sampling and finite size effects. In FigJ5]a 
and b (7 = 0.25 and 0.125) we have plotted the variation 
of JijS 2 as a function of the distance between impurities, 
for different system sizes (N = L 3 where L varies from 4 
to 16). We observe that the couplings are larger for the 
smallest systems and that finite size effects are especially 
strong at distances near the typical and relevant one, d. 
Thus, we already expect the Curie temperature to be 
strongly size dependent. In Fig|5]thc average Curie tem- 
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FIG. 8: (Color online.) Magnetic couplings Jij(r)S 2 in unit of 
t as a function of the the distance between impurities r/a.The 
magnetic impurity concentration is fixed x — 0.065, JS=5 t 
and (a) 7 = 0.125 and (b) 7 = 0.25 . The calculations are 
done with clusters of size N = L 3 where L= 4,8,10,12,14,16. 



perature as a function of JS/W is shown for various sizes. 
The hole density (per carrier) is set to 7 = 0.125. Note 
that each averaged T c is obtained using the associated 
couplings, e.g calculated for the corresponding system. 
Note also that the average was properly performed, we 
have used a sufficiently large number of disorder configu- 
rations (few thousands for the smallest systems and few 
hundreds for the largest). We analyze now the impor- 
tance of the statistical sampling. For a given value of 
JS, we observe beyond N = 14 3 an insignificant varia- 
tion of the averaged Curie temperature: beyond this size 
the thermodynamic limit is properly described. Because 
within standard full Monte Carlo method the systems are 
usually restricted to relatively small size a 32 i 33 ' 34 i 35 ' 36 of 
the order of N = 4 3 , let us do the comparison between the 
calculations done with the smallest and largest systems. 
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FIG. 9: (Color online.) Averaged Curie temperature (in units 
of the bandwidth W) as a function of the parameter JS/W. 
The magnetic impurity concentration is fixed x — 0.065. The 
hole concentration per impurity is 7 = 0.125. The calcu- 
lations are done with clusters of size N = L 3 where L= 
5,6,12,16,20. A systematic average over few thousands of con- 
figurations have been done for the smaller systems and few 
hundreds for the largest. 
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FIG. 10: (Color online.) Distribution of the calculated Curie 
temperatures for different systems sizes N=L 3 where L varies 
from 4 to 20. The parameters are x = 0.065, V — 0, 7 = 0.25 
and JS = 0.5 W 



For instance, we consider JS = 0.7 W. We see that the 
Curie temperature for L = 5 is Tc = 10.5 10~ 4 W whilst 
in the thermodynamic limit T c = 2.0 10" 4 W (500 % 
difference!). 

We now show that the sampling over disorder is also 
crucial and should be done properly, especially in di- 
lute systems. Note for instance that, Franceschetti et 
ali2£, in their ab initio based study, have shown impor- 
tant fluctuations of the Curie temperature from sample 
to sample. In Fig[TTJ] we have plotted the distribution 
of the Curie temperatures obtained using different sys- 
tem sizes. To facilitate the discussion, we have plot- 
ted on FigfTT] the averaged Curie temperatures, = 
j+oc Tc p L (j> c )dT c , and the w idths at mid height of the 
distributions (ATc/W) as a function of 1/N. For the 
smallest system (N = 4 3 ) one observes a very important 
spreading of the critical temperatures distribution: they 
can vary from at least one order of magnitude! From 
Fig[TTJ for this case, ATq/W is rather close to the cor- 
responding average value. For these small systems, one 
understands easily that 10 configurations of disorder are 
definitely insufficient to provide a reliable average. In 
other words, the calculated averaged Curie temperature 
over a small number of configurations (few tens) could 
easily lead to critical temperatures 10-20 times larger 
than that calculated properly (averaged over few thou- 
sands of configurations). This implies, for example, that 
the Curie temperatures shown in Figf9] could be easily 
50 times larger than the calculated ones in the thermo- 
dynamic limit if the average would have been taken over 
only few configurations. 
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FIG. 11: (Color online.) Width (AT C ) and average (T c ) 
of the Curie temperatures distribution calculated for various 
simple cubic lattices of size N=L 3 as a function of 1/N. The 
parameters are the same as those of Fig llOl 



As the system size increases, one gradually observes a 
decrease in ATc, which then rapidly tends to zero (ther- 
modynamic limit), whilst Tq tends towards a constant 
value (see FigfTT]). One can already consider that the 
thermodynamic limit value of Tc is reached for L > 16. 
In the light of the results presented in this manuscript, a 
criterion to get reasonable results would consist in con- 
sidering system sizes where L is at least five to six times 
larger than the typical distance between impurities. This 
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criterion is only rough, in the case of long range magnetic 
couplings, such as RKKY exchanges, the finite size effects 
can be even more drastic. In addition, only ten configu- 
rations of disorder for small systems are definitely insuffi- 
cient to get the average value of Tq in a reliable manner. 
It is obvious that these numerical requirements are very 
difficult to fulfill within standard Monte Carlo calcula- 
tions but they are definitely essential. In addition, in the 
presence of inhomogeneities, temperatures are also ex- 
pected to be even more sensitive to both finite size effects 
and statistical sampling 39 . Moreover, the way of extract- 
ing the Curie temperature is essential. Within SC-LRPA, 
the problem does not arise. The critical temperatures are 
directly given by a semi-analytical equation solved in a 
self-consistent way, no extrapolation from the magneti- 
zation curve is used. Concerning the Monte Carlo sim- 
ulations of reference^, the method used to extract Tc 
from the magnetization curve is not very accurate (see 
figure 4. a and 4.b of reL22). Potential errors related to 
this method are added to the uncertainties coming from 
both finite size effects and from the insufficiency of sta- 
tistical sampling. A more accurate way to obtain the 
Curie temperatures in MC simulations is based on the 
method of Binder cumulants 3 -^, but this method repre- 
sents an additional cost in terms of computing time since 
it needs a finite size effect analysis. Unfortunately, in 
spite of the small system sizes and the small number of 
configurations of disorder considered in r eiM, the Monte 
Carlo simulations already needs a huge resource in term 
of CPU time and memory. 

In a recent paper based on Monte Carlo simulations, 
Yildirim et al. have studied the diluted Kondo model ap- 
plied to the case of Gai-^Mn^As by including the realis- 
tic band structure of the host material 3 ^ (see also ref4^ 
for a comment). They have calculated the Curie temper- 
ature as a function of the local coupling J for systems con- 
taining x = 8.5% of magnetic impurities and ph = 0.75 
hole per impurity. We underline here that their numeri- 
cal Monte Carlo approach suffers from the same numeri- 
cal shortcomings as those previously discussed. The sys- 
tem considered in this letter contain typically 20 localized 
spins only and the average are done using just five con- 
figurations of disorder. Additionally, the value of J con- 
sidered in refill, corresponds to the perturbative regime. 
Thus, on the basis of the present study, we argue that 
within this limit, the couplings should exhibit RKKY os- 
cillations which should lead in the thermodynamic limit 
to small Curie temperatures or eventually to no ferromag- 
netic phased. Unfortunately, and as previously noticed, 
the smallness of the cluster considered in the MC cal- 
culations hide the effects of the asymptotic RKKY tail, 
and thus leading to finite and large Curie temperatures. 
We argue that by improving the statistics and increas- 



ing the size of the systems the Curie temperature should 
vanish when the frustration will become effective. Fur- 
thermore, we argue that the large MC critical temper- 
atures found in the large JS regime (see Fig. I of refill) 
are in fact a numerical artefact due again to the insuffi- 
cient statistical sampling and finite size effects. Indeed, in 
this regime double exchange mechanism dominates and 
thus very small or vanishing Curie temperatures are ex- 
pected, see the corresponding densities of states in Fig. 4 
of2i where the IB is separated from the valence band. 



V. CONCLUSION 

To conclude, in this work we have shown that the study 
of the diluted magnetic systems requires a rigorous nu- 
merical treatment. We have shed light on the origins of 
the disagreements between MC and TSA for the study of 
diluted magnetic systems. In contrast to the non dilute 
systems, the finite size effects as well as the importance 
of the statistical sampling appear to be crucial. We have 
shown that available Monte Carlo simulations for the di- 
luted model Hamiltonian ((T|) suffer from severe numerical 
insufficiencies. Although Monte Carlo simulations arc in 
principle exact, but because the calculations are unfor- 
tunately restricted to relatively small systems and weak 
statistical sampling, the obtained Curie temperatures are 
often largely overestimated. In addition, limiting regimes 
(strong and weak couplings) are not properly described. 
It would be of great interest to check both the importance 
of statistical sampling and the finite size effects within 
new large scale Monte Carlo studies. Among MC simu- 
lations which allow to reach relatively large system sizes, 
one can quote the Hybrid Monte Carlo method (HMC)^, 
the Polynomial Expansion Monte Carlo (PEMC)^ or 
its fastest counterpart, namely the Truncated Polyno- 
mial Expansion Monte Carlo (TPEMC) methodi^^iil. 
In Particular TPEMC handles higher systems up to 20 3 
sites within a reasonable CPU time (see table II of ref4£). 
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